%% Solution series 4 - Analysis of the design for the Menchutkin's experiments
%
% Table of Contents
% Design analysis
%     Loarding data
%     Data visualization
%     Computing the dispersion matrices for each model
%         Standardizing the eaasy matrix
%         A zone without data
%         Scatter matrices of the different models
%         VIF for the three models
%         Variance function
%         Analysis
% Optimization
%     The 15 best experiences from a D-optimal perspective
%         Dispersion matrix for 15 runs
%         Computing the VIF
%     The 10 best experiences from a D-optimal perspective
%         Dispersion matrix
%         Computing the VIF
%     The 5 best experiences from a D-optimal perspective
%         Dispersion matrix
%         Computing the VIF
% Degradation of the dispersion matrix
%         Evolution of the diagonal elements for the linear model
%         Evolution of the diagonal elements for the linear model with interactions
% Variance function
% Dispersion matrix and VIF
%
%% Design analysis
%
%% Loarding data
%
% The readtable() routine is used to load data from a file. The routine allows different types of files (.txt, .dat, .xlsx,...). It is also possible to use the Import Data utility found in the Home menu.
%

data=readtable('DataMenshutkin.xlsx','sheet','Menchutkin');
data.Amine=categorical(data.Amine);
disp(data)

%
%% Data visualization
%

figure
plot(data.sigmaF,data.sigmaR,'ok',...
    'MarkerFaceColor','r',...
    'MarkerSize',10)
set(gca,'FontName','Times New Roman',...
    'FontSize',14)
axis([0 1 -1 0.25])
title('Points de mesure',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('\sigma_F',...
    'FontName','Times New Roman',...
    'FontSize',18)
ylabel('\sigma_R',...
    'FontName','Times New Roman',...
    'FontSize',18)
grid on
text(data.sigmaF+.01,data.sigmaR-.03,data.Amine,...
    'FontName','Times New Roman',...
    'FontSize',14)

%
% We observe that there is a zone at the bottom right with no available amine, therefore no measurement point. There are also amines whose characteristics are very similar, so this will correspond to very close measurement points which will provide information little differentiated and therefore of less importance at first.
%
%% Computing the dispersion matrices for each model
%
% Even though experts in the field understand that there is a relationship between the dynamics of the reaction and the resonance and induction coefficients, there is no theoretical model available. We are therefore interested in empirical models. A priori a second degree model would be likely to indicate an optimum. However, the measurement points being imposed by the existing amines, it makes sense to also consider the simpler models. This is why the analysis will consider the three standard empirical models:
% linear R(\sigma_R,\sigma_F)=a_o+a_R \sigma_R+ a_F \sigma_F
% linear with interactions R(\sigma_R,\sigma_F)=a_o+a_R \sigma_R+ a_F \sigma_F + a_{RF} \sigma_R \sigma_F)
% quadraticR(\sigma_R,\sigma_F)=a_o+a_R \sigma_R+ a_F \sigma_F + a_{RF} \sigma_R \sigma_F) + a_{RR} \sigma_R^2 +a_{FF} \sigma_F^2
%% Standardizing the eaasy matrix

E=data{:,2:3}; % selecting the essay matrix
E=rescale(E,-1,1,"InputMin",min(E),"InputMax",max(E)); % standardized essay matrix

%
% Matlab also offers an alternative with the scatter() function with a certain number of graphic options that do not exist with plot() (see the function's help).

scatter(E(:,1),E(:,2),'red',"filled")
xlabel('x_F',...
    'FontName','Times New Roman',...
    'FontSize',18)
ylabel('x_R',...
    'FontName','Times New Roman',...
    'FontSize',18)
grid on

%%
%% A zone without data
% we can draw a line that serves as a boundary between the standardized coordinates (-0.6,-1) and (1.08)

hold on
    plot([-0.6;1],[-1;0.8],'-b')

% The coefficients \alpha of this line are obtained by computing  \alpha=A^{-1}R , if A is the model matrix

A=[1 -0.6;1 1];
disp(A)

% and R the corresponding values of x_2

R=[-1;0.8];
disp(R)
alpha=A\R

% So the line 45x_1-40x_2-13=0 is the border that we are looking for. The point (0,0) is in the ahthorized zone and correspondy to  a negative value -13. So the non-othautized zone is described by the inequality  45x_1-40x_2-13>0

patch([-0.6;1;1;-0.6],[-1;0.8;-1;-1],[.9 .9 .9])
hold off

%
%%
%% Scatter matrices of the different models
% As the calculation of the useful elements for the analysis of each model is repetitive, it is interesting to develop a function which performs these calculations automatically.

[D_lin,VIF_lin]=analyse(E,'linear');
[D_int,VIF_int]=analyse(E,'interactions');
[D_quad,VIF_quad]=analyse(E,'quadratic');

var_lin_19=diag(D_lin);
var_int_19=diag(D_int);

%

disp('')
disp('Modèle linéaire:')
fprintf('%5.2f\t%5.2f\t%5.2f\t\n',D_lin')
disp(' ')
disp('Modèle avec interaction:')
disp(' ')
fprintf('%5.2f\t%5.2f\t%5.2f\t%5.2f\t\n',D_int')
disp(' ')
disp('Modèle quadratique:')
disp(' ')
fprintf('%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t%5.2f\t\n',D_quad')

%
%% VIF for the three models

VIF=table;
VIF.Quadratique=VIF_quad;
VIF.Interaction(1:4,1)=VIF_int;
VIF.Lineaire(1:3,1)=VIF_lin;

figure
bar(VIF{:,1:3})
xlabel('Coefficients',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('VIF',...
    'FontName','Times New Roman',...
    'FontSize',14)
set(gca,'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_{12}' 'a_{11}' 'a_{22}'})
% grid on
hold on
plot([0;7],[6;6],':b',[0;7],[10;10],':r','LineWidth',2)
legend('Quadratique','Interaction','linéaire','Warning','Limite')
hold off

%
%% Variance function
% To simplify the calculation, which is repetitive, a function has been developed to plot the variance function for each case.

comp_var(D_lin,'linear');
comp_var(D_int,'interactions');
comp_var(D_quad,'quadratic');

%
%% Analysis
% The values ��of the diagonal elements of the dispersion matrices are acceptable for the linear model, limits for the model with interaction and frankly alarming for the quadratic model.
% VIF values ��follow the same trend
% The variance functions clearly show that the linear model is acceptable for the entire domain while the functions for the other two models clearly show significant problems in the small and large area which corresponds to the area with no measurement point.
% In this case, it is therefore on the dispersion matrix and the VIFs that the analysis must be based. For the variance function, the fact that the model becomes very imprecise in the zone where it is extrapolated is not a criterion in itself because the model will not be used in this zone. The function should be cut using nan (not a number). Here is the example with the interaction model:

syms x1 x2
f=[1; x1; x2; x1*x2];
fvar=f.'*D_int*f;
[xx1,xx2]=meshgrid(-1:.1:1);
fvar_num=double(subs(fvar,{x1,x2},{xx1,xx2}));
k=find(45*xx1-40*xx2-13>0);
fvar_num(k)=nan;
surf(xx1,xx2,fvar_num)

view([-73.90 65.81])
title(['Fonction de variance: modèle avec interaction'],...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('x_1',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('x_2',...
    'FontName','Times New Roman',...
    'FontSize',14)
zlabel('var/s^2',...
    'FontName','Times New Roman',...
    'FontSize',14)

xlim([-1.00 1.00])
ylim([-1.00 1.00])
zlim([0.00 2.35])
view([-75.77 37.09])

%
%
%% Optimization
%% The 15 best experiences from a D-optimal perspective
% The candexch() function allows to select the best D-optimal experiments.

list_15=candexch(x2fx(E,'interaction'),15);
table_15=table(sort(list_15),'VariableNames',{'Expérience'});
table_15.Amine=data.Amine(sort(list_15));
table_15.sigmaF=data.sigmaF(sort(list_15));
table_15.sigmaR=data.sigmaR(sort(list_15));

disp(table_15)

%
% We can represent the measurement points on the experimental space and observe that the algorithm chooses the most distant points and prefers to duplicate good experiments than to choose other points. The choice of points will depend on the model. Try other models!
%

figure
plot(data.sigmaF(list_15,:),data.sigmaR(list_15,:),'ok',...
    'MarkerFaceColor','r',...
    'MarkerSize',10)
set(gca,'FontName','Times New Roman',...
    'FontSize',14)

title('15 meilleurs points de mesure pour un modèle avec interaction',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('\sigma_F',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('\sigma_R',...
    'FontName','Times New Roman',...
    'FontSize',14)
grid on
text(data.sigmaF(list_15,:)+.03,data.sigmaR(list_15,:)-.02,data.Amine(list_15,:),...
    'FontName','Times New Roman',...
    'FontSize',14)

%

% matrice d'expériences standardisée
E15=data{list_15,2:3}; 
E15=rescale(E15,-1,1,"InputMin",min(E15),"InputMax",max(E15));

[D_lin,VIF_lin]=analyse(E15,'linear');
[D_int,VIF_int]=analyse(E15,'interactions');
var_lin_15=diag(D_lin); % met en réserve les éléments diagonaux de D
var_int_15=diag(D_int);

%
%% Dispersion matrix for 15 runs

disp(' ')
disp('Modèle linéaire:')
fprintf('%5.2f\t%5.2f\t%5.2f\t\n',D_lin')
disp(' ')
disp('Modèle avec interaction:')
disp(' ')
fprintf('%5.2f\t%5.2f\t%5.2f\t%5.2f\t\n',D_int')

%
%% Computing the VIF

clear VIF
VIF=table;
VIF.Interaction(1:4,1)=VIF_int;
VIF.Lineaire(1:3,1)=VIF_lin;

figure
bar(VIF{:,1:2})
title('Avec 15 essais',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('Coefficients',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('VIF',...
    'FontName','Times New Roman',...
    'FontSize',14)
set(gca,'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_{12}'})
% grid on
hold on
plot([0;5],[6;6],':b',[0;5],[10;10],':r','LineWidth',2)
legend('Interaction','linéaire','Warning','Limite')
hold off

%
%% The 10 best experiences from a D-optimal perspective

list_10=candexch(x2fx(E,'interaction'),10);
table_10=table(sort(list_10),'VariableNames',{'Expérience'});
table_10.Amine=data.Amine(sort(list_10));
table_10.sigmaF=data.sigmaF(sort(list_10));
table_10.sigmaR=data.sigmaR(sort(list_10));

disp(table_10)

figure
plot(data.sigmaF(list_10,:),data.sigmaR(list_10,:),'ok',...
    'MarkerFaceColor','r',...
    'MarkerSize',10)
set(gca,'FontName','Times New Roman',...
    'FontSize',14)

title('10 meilleurs points de mesurepour un modèle avec interaction',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('\sigma_F',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('\sigma_R',...
    'FontName','Times New Roman',...
    'FontSize',14)
grid on
text(data.sigmaF(list_10,:)+.03,data.sigmaR(list_10,:)-.02,data.Amine(list_10,:),...
    'FontName','Times New Roman',...
    'FontSize',14)

% matrice d'expériences standardisée
E10=data{list_10,2:3}; 
E10=rescale(E10,-1,1,"InputMin",min(E10),"InputMax",max(E10));


[D_lin,VIF_lin]=analyse(E10,'linear');
[D_int,VIF_int]=analyse(E10,'interactions');

var_lin_10=diag(D_lin); % met en réserve les éléments diagonaux de D
var_int_10=diag(D_int);

%% Dispersion matrix

disp(' ')
disp('Modèle linéaire:')
fprintf('%5.2f\t%5.2f\t%5.2f\t\n',D_lin')
disp(' ')
disp('Modèle avec interaction:')
disp(' ')
fprintf('%5.2f\t%5.2f\t%5.2f\t%5.2f\t\n',D_int')

%% Computing the VIF

VIF=table;
VIF.Interaction(1:4,1)=VIF_int;
VIF.Lineaire(1:3,1)=VIF_lin;

figure
bar(VIF{:,1:2})
title('Avec 10 essais',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('Coefficients',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('VIF',...
    'FontName','Times New Roman',...
    'FontSize',14)
set(gca,'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_{12}'})
% grid on
hold on
plot([0;5],[6;6],':b',[0;5],[10;10],':r','LineWidth',2)
legend('Interaction','linéaire','Warning','Limite')
hold off

%
%% The 5 best experiences from a D-optimal perspective

list_5=candexch(x2fx(E,'interaction'),5);
table_5=table(sort(list_5),'VariableNames',{'Expérience'});
table_5.Amine=data.Amine(sort(list_5));
table_5.sigmaF=data.sigmaF(sort(list_5));
table_5.sigmaR=data.sigmaR(sort(list_5));
disp(table_5)

figure
plot(data.sigmaF(list_5,:),data.sigmaR(list_5,:),'ok',...
    'MarkerFaceColor','r',...
    'MarkerSize',10)
set(gca,'FontName','Times New Roman',...
    'FontSize',14)

title('5 meilleurs points de mesure',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('\sigma_F',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('\sigma_R',...
    'FontName','Times New Roman',...
    'FontSize',14)
grid on
text(data.sigmaF(list_5,:)+.03,...
    data.sigmaR(list_5,:)-.02,...
    data.Amine(list_5,:),...
    'FontName','Times New Roman',...
    'FontSize',14)

% matrice d'expériences standardisée
E5=data{list_5,2:3};
E5=rescale(E5,-1,1,"InputMin",min(E5),"InputMax",max(E5));
[D_lin,VIF_lin]=analyse(E5,'linear');
[D_int,VIF_int]=analyse(E5,'interactions');

var_lin_5=diag(D_lin); % met en réserve les éléments diagonaux de D
var_int_5=diag(D_int);

%
%% Dispersion matrix

disp(' ')
disp('Modèle linéaire:')
fprintf('%5.2f\t%5.2f\t%5.2f\t\n',D_lin')
disp(' ')
disp('Modèle avec interaction:')
disp(' ')
fprintf('%5.2f\t%5.2f\t%5.2f\t%5.2f\t\n',D_int')

%% Computing the VIF

VIF=table;
VIF.Interaction(1:4,1)=VIF_int;
VIF.Lineaire(1:3,1)=VIF_lin;

figure
bar(VIF{:,1:2})
title('Avec 5 essais',...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('Coefficients',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('VIF',...
    'FontName','Times New Roman',...
    'FontSize',14)
set(gca,'XTickLabel',{'a_o' 'a_1' 'a_2' 'a_{12}'})
% grid on
hold on
plot([0;5],[6;6],':b',[0;5],[10;10],':r','LineWidth',2)
legend('Interaction','linéaire','Warning','Limite')
hold off

%%
%% Degradation of the dispersion matrix
%
%% Evolution of the diagonal elements for the linear model

var_lin=table(var_lin_19,var_lin_15,var_lin_10,var_lin_5,...
    'RowName',{'a_o' 'a_1' 'a_2' } );
disp(var_lin)

% The diagonal elements of the dispersion matrix remain acceptable even with only 5 experiments with the variance of the coefficients which is then multiplied approximately three times.
%%
%% Evolution of the diagonal elements for the linear model with interactions

var_int=table(var_int_19,var_int_15,var_int_10,var_int_5,...
    'RowName',{'a_o' 'a_1' 'a_2' 'a_12'}  );
disp(var_int)

% The diagonal elements of the dispersion matrix quickly become too large if the number of experiments is reduced. Limiting the number of experiments to 10 already produces coefficients with variance equal to the experimental variance.
%
%
%% Variance function

function comp_var(D,spec)
    % * Définition des variables symboliques pour les 3 variables
    syms x1 x2 ;

    if strcmp(spec,'linear')
        f=[1; x1; x2];
    elseif strcmp(spec,'interactions')
        f=[1; x1; x2; x1*x2];
    elseif strcmp(spec,'quadratic')
        f=[1; x1; x2; x1*x2; x1^2; x2^2];
    else
         error('Erreur dans la fonction de variance: Le modèle n''a pas été reconnu') 
    end

    % Calcul de la fonction de variance
    fvar=f.'*D*f;
    
    figure
    
fsurf(fvar, [-1 1 -1 1],'MeshDensity',16)
title(['Fonction de variance: ' spec],...
    'FontName','Times New Roman',...
    'FontSize',16)
xlabel('x_1',...
    'FontName','Times New Roman',...
    'FontSize',14)
ylabel('x_2',...
    'FontName','Times New Roman',...
    'FontSize',14)
zlabel('var/s^2',...
    'FontName','Times New Roman',...
    'FontSize',14)
     
    
end

%%
%% Dispersion matrix and VIF

function [D,VIF]=analyse(E,spec)
X=x2fx(E,spec);
D=inv(X'*X);
C=corrcov(D);
VIF=diag(inv(C));

end